Derived quantities we would like to compute uncertainty for:
$C_P = \frac{T \omega}{\frac{1}{2} \rho A U_\infty^3}$
$C_D = \frac{F_L + F_R}{\frac{1}{2} \rho A U_\infty^2}$
$\lambda = \frac{\omega R}{U_\infty}$
We will use the methods defined in Section 2-5 of Coleman and Steele (2009), where uncertainty is defined as
$$ u_c^2 = s_X^2 + \sum_{k=1}^{M} b_k^2, $$where $u$ is the uncertainty, $s_X$ is the sample standard deviation, and $b_k$ is a source of systematic error, which can be obtained from an instrument's calibration/datasheet. For example, if a datasheet claims an instrument to be accurate within $\pm 1.0 \text{ N}$, this would be equal to $2b$, and therefore $b = 0.5 \text{ N}$.
Note that for our quantities we are looking for uncertainties of the means, for which the standard deviation is $$ s_{\bar{X}} = \frac{1}{\sqrt{N}} s_X, $$ where $N$ is the number of independent samples, which in our case we choose to be turbine revolutions.
| Quantity | Nominal accuracy | $b$ |
|---|---|---|
| Torque | $\pm 0.5 \text{ Nm}$ | $0.25 \text{ Nm} $ |
| Turbine angle | $2\pi/10^{5}/2 = \pm 3 \times 10^{-5} \text{ rad}$ | $1.5 \times 10^{-5} \text{ rad}$ |
| Carriage position | $\pm 5 \times 10^{-6} \text{ m}$ | $2.5 \times 10^{-6} \text{ m}$ |
| Drag force (one side) | $\pm 0.06 \text{ lbf} = 0.28 \text{ N} $ | $0.14 \text{ N} $ |
In [3]:
# Get things all setup
from __future__ import division, print_function
import os
import numpy as np
import scipy.stats
if os.path.split(os.getcwd())[-1] == "IPython notebooks":
print("Setting working directory to experiment root directory")
os.chdir("../../")
from pyrvatrd.processing import *
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline
plt.style.use("ggplot")
b_torque = 0.25
b_angle = 1.5e-5
b_car_pos = 2.5e-6
b_force = 0.14
In [4]:
def calc_uncertainty(quantity, b):
"""Calculate uncertainty of a mean quantity."""
n = len(quantity)
return np.sqrt((quantity.std()/np.sqrt(n))**2 + b**2)
In [6]:
run = Run("Wake-1.0", 20)
torque = run.torque
unc_torque = calc_uncertainty(torque, b_torque)
print("Torque uncertainty (using all samples) = {:0.1e} Nm".format(unc_torque))
torque = run.torque_per_rev
unc_torque = calc_uncertainty(torque, b_torque)
print("Torque uncertainty (averaging per rev) = {:0.1e} Nm".format(unc_torque))
From Coleman and Steele (2009) equation 3.15, the uncertainty of the power coefficient is
$$ u_{C_P}^2 = s_{C_P}^2 + b_{C_P}^2, $$where
$$ b_{C_P}^2 = \sum_{i=1}^J \left( \frac{\partial C_P}{\partial X_i} \right)^2 b_{X_i}^2. $$For our case, we consider $X = [T, \omega, U_\infty ] $. The partial derivatives are then
$$ \frac{\partial C_P}{\partial T} = \frac{\omega}{\frac{1}{2} \rho A U_\infty^3}, $$$$ \frac{\partial C_P}{\partial \omega} = \frac{T}{\frac{1}{2} \rho A U_\infty^3}, $$$$ \frac{\partial C_P}{\partial U_\infty} = \frac{-3 T \omega}{\frac{1}{2} \rho A U_\infty^4}. $$
In [7]:
rho = 1000.0
area = 1.0
omega = run.omega.mean()
torque = run.torque.mean()
u_infty = run.tow_speed.mean()
const = 0.5*rho*area
b_cp = np.sqrt((omega/(const*u_infty**3))**2*b_torque**2 + \
(torque/(const*u_infty**3))**2*b_angle**2 + \
(-3*torque*omega/(const*u_infty**4))**2*b_car_pos**2)
unc_cp = calc_uncertainty(run.cp_per_rev, b_cp)
print("Uncertainty of C_P = {:.3f}".format(unc_cp))
Drag coefficient is calculated using
$$ C_D = \frac{F_L + F_R}{\frac{1}{2} \rho A U_\infty^2} . $$The measured variables used to calculate $C_D$ are $[F_L, F_R, U_\infty]$, therefore the necessary partial derivatives are
$$ \frac{\partial C_D}{\partial F_L} = \frac{1}{\frac{1}{2} \rho A U_\infty^2} ,$$$$ \frac{\partial C_D}{\partial F_R} = \frac{1}{\frac{1}{2} \rho A U_\infty^2} ,$$$$ \frac{\partial C_D}{\partial U_\infty} = \frac{-2(F_L + F_R)}{\frac{1}{2} \rho A U_\infty^3} .$$
In [8]:
drag = run.drag.mean()
u_infty = run.tow_speed.mean()
const = 0.5*rho*area
b_cd = np.sqrt((1/(const*u_infty**2))**2*b_force**2 + \
(1/(const*u_infty**2))**2*b_force**2 +
(-2*drag/(const*u_infty**3))**2*b_car_pos**2)
unc_cd = calc_uncertainty(run.cd_per_rev, b_cd)
print("Uncertainty of C_D = {:.3f}".format(unc_cd))
The necessary partial derivatives are
$$ \frac{\partial \lambda}{\partial \omega} = \frac{R}{U_\infty} , $$$$ \frac{\partial \lambda}{\partial U_\infty} = \frac{-\omega R}{U_\infty^2} . $$
In [9]:
omega = run.omega.mean()
u_infty = run.tow_speed.mean()
R = 0.5
b_tsr = np.sqrt((R/(u_infty))**2*b_angle**2 + \
(-omega*R/(u_infty**2))**2*b_car_pos**2)
unc_tsr = calc_uncertainty(run.tsr_per_rev, b_tsr)
print("Uncertainty of TSR = {:.3f}".format(unc_tsr))
The uncertainty calculated above is the combined standard uncertainty. If we want to apply a confidence level, we will obtain the expanded uncertainty
$$ U_{\%} = k_{\%}u_c ,$$where $k_{\%}$ is usually chosen from the Student-$t$ distribution for a given percent level of confidence. To choose a $t$-value, we need an estimate for degrees of freedom, which can be obtained from the Welch--Satterthwaite formula
$$ \nu_X = \frac{\left(s_X^2 + \sum_{k=1}^M b_k^2 \right)^2} {s_X^4/\nu_{s_X} + \sum_{k=1}^M b_k^4/\nu_{b_k}}, $$where $\nu_{s_X}$ is the number of degrees of freedom associated with $s_X$ and $\nu_{b_k}$ is the number of degrees of freedom associated with $b_k$. $\nu_{s_X}$ is simply $N-1$, where the ISO guide recommends
$$ \nu_{b_k} = \frac{1}{2} \left( \frac{\Delta b_k}{b_k} \right)^{-2}, $$where the quantity in parentheses is the relative uncertainty of $b_k$.
In [11]:
# Let's try to calculate nu_cp
s_cp = run.cp_per_rev.std()
nu_s_cp = len(run.cp_per_rev) - 1
# Already have b_cp from above
b_cp_rel_unc = 0.25 # A guess
nu_b_cp = 0.5*b_cp_rel_unc**(-2)
nu_cp = ((s_cp**2 + b_cp**2)**2)/(s_cp**4/nu_s_cp + b_cp**4/nu_b_cp)
t = scipy.stats.t.interval(alpha=0.95, df=nu_cp)[-1]
exp_unc_cp = t*unc_cp
print("Degrees of freedom nu_cp = {:.0f}".format(nu_cp))
print("Expanded uncertainty of C_P (95% conf) = {:.2f}".format(exp_unc_cp))
# Let's try to calculate nu_cd
s_cd = run.std_cd_per_rev
nu_s_cd = len(run.cd_per_rev) - 1
# Already have b_cd from above
b_cd_rel_unc = 0.25 # A guess
nu_b_cd = 0.5*b_cd_rel_unc**(-2)
nu_cd = ((s_cd**2 + b_cd**2)**2)/(s_cd**4/nu_s_cd + b_cd**4/nu_b_cd)
t = scipy.stats.t.interval(alpha=0.95, df=nu_cd)[-1]
exp_unc_cd = t*unc_cd
print("Degrees of freedom nu_cd = {:.0f}".format(nu_cd))
print("Expanded uncertainty of C_D (95% conf) = {:.2f}".format(exp_unc_cd))
# Let's try to calculate nu_tsr
s_tsr = run.std_tsr_per_rev
nu_s_tsr = len(run.tsr_per_rev) - 1
# Already have b_tsr from above
b_tsr_rel_unc = 0.25 # A guess
nu_b_tsr = 0.5*b_tsr_rel_unc**(-2)
nu_tsr = ((s_tsr**2 + b_tsr**2)**2)/(s_tsr**4/nu_s_tsr + b_tsr**4/nu_b_tsr)
t = scipy.stats.t.interval(alpha=0.95, df=nu_tsr)[-1]
exp_unc_tsr = t*unc_tsr
print("Degrees of freedom nu_tsr = {:.0f}".format(nu_tsr))
print("Expanded uncertainty of TSR (95% conf) = {:.3f}".format(exp_unc_tsr))